********************************
** Read in file, add controls **

use 	"./output/citiesall_patents_ageyoungold.dta", clear

merge 1:1 gen_id f_myr using "./output/citiesall_patents_migration.dta"
drop _merge

browse gen_id city f_myr master_id pat_wtd_invadj pat_overeq40 pat_under40 pat_wtd_inv totpat_stayer*

merge 1:1 gen_id f_myr using "./output/cities_expandedpatents_breakthrough"
drop _merge


** npi X month-of-year FEs **
egen	idxm = group(gen_id f_m)

** census region X year FEs **
egen	regyr = group(f_yr region)

keep if tin(1915m1, 1925m12)
*=======================*
*= Treatment variables  *
*=======================*

tab 	days_npis if f_myr==tm(1918m9)

gen		daysnpis = days_npis/30 
gen		lnnpis = ln(days_npis)
gen		longnpis = (days_npis>=90) if !mi(days_npis) // Natural break in data, about 61% have short
gen		ldeaths = ln(mag_maxexcessdeaths)

foreach v of varlist daysnpis longnpis lnnpis ldeaths {
	gen 	post_X_`v' 		= `v'*(f_myr>=tm(1918m9))

	gen 	during_X_`v' 	= `v'*( (f_myr>=tm(1918m9)) & (f_myr<=tm(1919m3)) )
	gen 	after_X_`v' 	= `v'*(f_myr>=tm(1919m4))
	gen 	pre_X_`v' 		= `v'*(f_myr<=tm(1917m8))
	
	lab var post_X_`v'		"Post Pandemic $\times$ NPI Length"
	
	lab var pre_X_`v'		"Before Pandemic $\times$ NPI Length"
	lab var during_X_`v' 	"During Pandemic $\times$ NPI Length"
	lab var after_X_`v' 	"After Pandemic $\times$ NPI Length"
}

drop 	daysnpis 

sort gen_id f_myr 

gen totpat_movers = pat_wtd_inv - totpat_stayers
gen totpat_internat = max(pat_wtd_inv - totpat_usborn, 0) // Only two cells with slightly negative values


*=======================*
*= Table 2, Panel C Regression Models =*
*=======================*

est clear

** Combination of  Models

foreach v in longnpis {
	
	local treat post_X_`v'

	eststo: ppmlhdfe pat_overeq40 `treat' if tin(1916m1, 1920m12), exp(pop_int) a(f_myr idxm regyr) cluster(gen_id)
		sum		pat_overeq40 if tin(1916m1, 1920m12)
		estadd scalar outmean = `r(mean)'
	eststo: ppmlhdfe pat_under40 `treat' if tin(1916m1, 1920m12), exp(pop_int) a(f_myr idxm regyr) cluster(gen_id)
		sum		pat_under40 if tin(1916m1, 1920m12)
		estadd scalar outmean = `r(mean)'
		
	eststo: ppmlhdfe totpat_stayers `treat' if tin(1916m1, 1920m12), exp(pop_int) a(f_myr idxm regyr) cluster(gen_id)
		sum		totpat_stayers if tin(1916m1, 1920m12)
		estadd scalar outmean = `r(mean)'	
	eststo: ppmlhdfe totpat_movers `treat' if tin(1916m1, 1920m12), exp(pop_int) a(f_myr idxm regyr) cluster(gen_id)
		sum		totpat_movers if tin(1916m1, 1920m12)
		estadd scalar outmean = `r(mean)'	
		
	eststo: ppmlhdfe totpat_usborn `treat' if tin(1916m1, 1920m12), exp(pop_int) a(f_myr idxm regyr) cluster(gen_id)
		sum		totpat_usborn if tin(1916m1, 1920m12)
		estadd scalar outmean = `r(mean)'	
	eststo: ppmlhdfe totpat_internat `treat' if tin(1916m1, 1920m12), exp(pop_int) a(f_myr idxm regyr) cluster(gen_id)
		sum		totpat_internat if tin(1916m1, 1920m12)
		estadd scalar outmean = `r(mean)'	
		
	eststo: ppmlhdfe break_p50_rrfsim05_wtd_inv `treat' if tin(1916m1, 1920m12), exp(pop_int) a(f_myr idxm regyr) cluster(gen_id)
		sum		break_p50_rrfsim05_wtd_inv if tin(1916m1, 1920m12)
		estadd scalar outmean = `r(mean)'		
	eststo: ppmlhdfe break_p50_rrfsim010_wtd_inv `treat' if tin(1916m1, 1920m12), exp(pop_int) a(f_myr idxm regyr) cluster(gen_id)
		sum		break_p50_rrfsim010_wtd_inv if tin(1916m1, 1920m12)
		estadd scalar outmean = `r(mean)'	
	
	estout 	using "$RES/fortable/combohet_posteffect_`v'.tex", replace cells(b(star fmt(%9.3f)) se(par)) ///
		stats(outmean N, fmt(%9.2f %9.0g) labels("Mean of Dep.\ Variable" "\(N\)"))  starlevels(+ 0.10 * 0.05 ** 0.01 *** 0.001) ///
		style(tex) drop(_cons) eqlabels(none) label
	est 	clear
}
